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Abstract In recent years, binary quadratic programming (BQP) has been 
successively applied to solve several combinatorial optimization problems. We 
consider in this paper a study of using the BQP model to solve the minimum 
sum coloring problem (MSCP). For this purpose, we recast the MSCP with 
a quadratic model which is then solved via a recently proposed Path Relink- 
ing (PR) algorithm designed for the general BQP. Based on a set of MSCP 
benchmark instances, we investigate the performance of this solution approach 
compared with existing methods. 
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1 Introduction 

Given an undirected graph G — (V, E) with vertex set V and edge set E, a 
K-coloring of G is a function c : v H> c(v) that assigns to each vertex v G V a 
color c(v), where c(v) G {1, 2, . . . , K}. A K-coloring is considered legal if each 
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pair of vertices (u, v) connected by an edge (u, v) £ E receive different colors 
c(u) ^ c(v). The minimum sum coloring problem (MSCP) is to find a legal 
K-coloring c such that the total sum of colors over all the vertices Ylvev c ( v ) 
is minimized. The minimum value of this sum is called the chromatic sum of 
G and denoted by ^2(G). The number of colors related to the chromatic sum 
is called the strength of the Graph and denoted by s(G). 

The MSCP is NP-hard for general graphs [T5] and provides applications 
mainly including VLSI design, scheduling and resource allocation [51114] . Given 
the theoretical and practical significance of the MSCP, effective approximation 
algorithms and polynomial algorithms have been presented for some special 
cases of graphs, such as trees, interval graphs and bipartite graphs [4TI51HU1 
IT2llT3"IIT4Tl2"l"ll2"2] . For the purpose of practical solving of the general MSCP, a 
variety of heuristics have been proposed in recent years, comprising a paral- 
lel genetic algorithm [TB], a greedy algorithm [T5], a tabu search algorithm 
[6], a hybrid local search algorithm [7], an independent set extraction based 
algorithm [25] and a local search algorithm [llj . 

On the other hand, binary quadratic programming (BQP) has emerged 
during the past decade as a unified model for a wide range of combinatorial 
optimization problems, such as set packing [2J, set partitioning [5D], general- 
ized independent set [17] . maximum edge weight clique [T] and maxcut |23) . 
A review concerning the additional applications and the reformulation proce- 
dures can be found in [IB]. This BQP approach has the advantage of directly 
applying an algorithm designed for BQP to solve other classes of problems 
rather than resorting to a specialized solution method. Moreover, this ap- 
proach proves to be competitive or even better than the special algorithms 
proposed for several problems. 

In this paper, we investigate for the first time the application of this BQP 
approach to solve the MSCP problem. We propose a binary quadratic formula- 
tion for the MSCP which is solved by our Path Relinking algorithm previously 
designed for the general BQP. To assess the performance of the proposed ap- 
proach, we present computational results on a set of 23 benchmark instances 
from the literature and contrast these results with those of several reference 
algorithms specifically dedicated to the MSCP. 

The rest of this paper is organized as follows. Section [5J illustrates how to 
transform the MSCP into the BQP formulation. Section[3]presents an overview 
of our Path Relinking algorithm for the general BQP. Section @] is dedicated 
to computational results and comparisons with other reference algorithms in 
the literature. The paper is concluded in Section [5] 



2 Transformation of the MSCP to the BQP model 

2.1 Linear model for the MSCP 

Given an undirected graph G = (V, E) with vertex set V (n = | V|) and edge 
set E. Let x u k be 1 if vertex u is assigned color k, and otherwise. The linear 
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programming model for the MSCP can be formulated as follows: 

n K 



Min f(x) = ^ ^ k ■ x uk 



u=l fc=l 
K 

subject to: cl. / j X u k = 1, uG{l,...,n} (-U 

fe=i 
c2. x ufe + x vk < 1, V(u, v) e E, k e {1, . . . , K} 
c3. x ufe e {0,1} 

2.2 Nonlinear BQP alternative 

The linear model of the MSCP can be recast into the form of the BQP ac- 
cording to the following steps: 

For the constraints cl., we represent these linear equations by a matrix 
Ax = b and incorporate the following penalty transformation |16j : 

#1: fl{x)=P(Ax-b) t {Ax-b) 

= P[x t {A t A)x - x*(A*&) - (b*A)x + b l b] 
= P[x t {A t A)x ~ a?*(A*6) - [b t A)x] + Ptfb 
= xD\x + c 



(2) 



69 For the constraints c2., we utilize the quadratic penalty function g(x) — 

70 PxukXyk to replace each inequality x u k + x v k < 1 in c2. and add them up as 

71 follows [16] : 

n n K 

#2: f2(x) = y^ y~^ y^ w uv x uk x vk 

u—l v—l,u^v k—1 
= XD2X 

(3) 

72 where w ut , = P if (u, v) £ E and otherwise. 

73 To construct the nonlinear BQP formulation h(x), we first inverse the 

74 minimum objective of the MSCP to be —f(x) in accordance with the general 

75 BQP model under a maximum objective, which becomes the first component 

76 of h(x). Then we add the penalty function fl(x) into h(x) such that fl(x) = 

77 if all the linear equations in cl. are satisfied and otherwise /l(x) is a penalty 

78 term with large negative values. In the same way, we add the penalty function 

79 f2(x) into h(x). Hence, the resulting BQP formulation for the MSCP can be 
so expressed as follows: 

BQP: Maxh(x) = -f(x) + fl(x)+f2(x) 

x + c 
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Once the optimal objective value for this BQP formulation is obtained, 
the minimum sum coloring value can be readily obtained by taking its inverse 
value. 

Further, a penalty scalar P is considered to be suitable as long as its 
absolute value \P\ is larger than half of the maximum color (\P\ > K/2). 
Consider that penalty functions should be negative under the case of a maximal 
objective, we select P = —500 for the benchmark instances experimented in 
this paper. The optimized solution x obtained by solving the nonlinear BQP 
formulation indicates that such selection ensures both fl(x) and f2(x) equal to 
0. In other words, each variable x u k with the assignment of 1 in the optimized 
solution x forms a feasible K-coloring in which vertex u gets the color k. 



92 2.3 An example of the transformation 

93 To illustrate the transformation from the MSCP to the BQP formulation, we 

94 consider the following graph with | V\ = 4 and expect to find a legal AT-coloring 
with K = 2. 




Fig. 1 A small graph sample 



Its linear formulation according to Equation ([l} is: 



Max f{x) = -x\\ - 2x\2 - #21 - 2x 2 2 - #31 - 2x 32 - x u - 2a; 4 2 
cl. 



c2. 



c3. 



#11 +#12 = 1 




#21 + #22 = 1 




#31 + #32 = 1 




#41 + #42 = 1 




#11 + #21 < 1 


#12 +2^22 < 1; 


#11 + #41 < 1 


#12 + #42 < 1; 


#21 +X 3 i < 1 


#22 +#32 < 1; 


#31 + #41 < 1 


#32 +#42 < 1; 


#11, #12, ■ • • ,# 


42 G {0,1} 
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Choosing the scalar penalty P = —5, we obtain the following BQP model: 
Max f{x) = xQx - 20 
where xQx is written as: 



Xn X12 



, . X41 X42 ) X 



/ 4 -5 -5 -5 \ 






-5 3 0-5 0-5 




( Xll \ 


-5 4-5-5 




X\2 


0-5-5 3 0-5 






0-5 4-5-5 


X 




0-5-5 3 0-5 




£41 


-5 0-5 4-5 


\X41J 


\ -5 0-5-53/ 







Solving this BQP model yields Xu — £22 = X31 = X42 = 1 (all other 
variables equal zero) and the optimal objective function value fix) = —6. 
Reversing this objective function value leads to the optimum (the minimum 
sum coloring) of 6 for this graph. 



104 
105 
106 
107 
108 
109 
110 
111 
112 
113 
114 
115 
116 
117 
118 
119 
120 
121 
122 
123 
124 
125 
126 



3 Path Relinking algorithm 

Path Relinking 9 is a general search strategy closely associated with tabu 
search and its underlying ideas share a significant intersection with the tabu 
search perspective [5], with applications in a variety of contexts where it has 
proved to be very effective in solving difficult problems. Our previous Path 
Relinking algorithm (PR) is directly utilized to solve the MSCP in its nonlin- 
ear BQP form as expressed in Eq. 4. The proposed Path Relinking algorithm 
is mainly composed of the following components: RefSet Initialization, Solu- 
tion Improvement, Relinking, Solution Selection, RefSet Updating and RefSet 
Rebuilding. Below we overview the general scheme of our Path Relinking al- 
gorithm. More details can be found in [2~4] . 

Firstly the RefSet Initialization method is used to create an initial set 
of elite solutions RefSet where each solution is obtained by applying a tabu 
search procedure to a randomly generated solution. 

Secondly, for each pair of solutions (x l , x 3 ) in RefSet we undertake the 
following operations: (1) apply a Relinking method to generate two paths 
(from x l to x 3 and from x 3 to x l ) by exploring trajectories (strictly confined 
to the neighborhood space) that connect high-quality solutions. The first and 
last solutions of the path are respectively called the initiating and guiding 
solutions. At each step of building the path from the initiating solution to the 
guiding solution, we randomly select a variable from the set of variables for 
which x l and x 3 ) have different values; (2) apply a Solution Selection method 
to select one solution from those generated on the path by reference both to its 
quality and to the hamming distance of this solution to the initiating solution 
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127 and the guiding solution. This selected solution is then submitted to the tabu 

128 search based Improvement Method; (3) apply a RefSet Updating method to 

129 decide if the newly improved solution is inserted in RefSet to replace the 
oo worst solution. The above procedure is called a round of our Path Relinking 
oi procedure. 

132 Finally, once a round of Path Relinking procedure is completed, a RefSet 

133 Rebuilding method is launched to rebuild RefSet that is used by the next round 

134 of the Path Relinking procedure. The Path Relinking algorithm terminates as 

135 soon as its stop condition (e.g. a fixed computing time) is satisfied. 



us 4 Experimental results 

137 4.1 Experimental protocols 

us To assess this BQP approach for the MSCP, we carry out experiments on a set 
139 of 23 graphtLj, which are the most used benchmark instances in the literature. 
wo Our experiments are conducted on a PC with Pentium 2.83GHz CPU and 
«i 8GB RAM. The time limit for a single run of our Path Relinking algorithm 
142 is set as follows: 1 hour for the first 16 instances in Table Q] 10 hours for 
dsjcl25.1, dsjcl25.5, dsjcl25.9, dsjc250.1 and dsjc250.5; 20 hours for dsjc250.9 
and dsjc500.1. Given the stochastic nature of our PR algorithm, each problem 
instance is independently solved 20 times. 

The tabu tenure til and the improvement cutoff \x |24] are two parameters in 
the tabu search based improvement method-a key component of the PR algo- 
rithm. According to preliminary experiments, we set til = max{40, A/100 + 
149 rand(50)} (where A denotes the number of variables in the resulting BQP 
model and rand(50) receives a random integer ranging from 1 to 50). In ad- 
dition, we set [i = 2A for the improvement of the initial solutions in RefSet 
and /i = 500 for the improvement of the solutions on the path, respectively. 
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4.2 Results of the BQP model for the MSCP 

Table Q] presents the computational statistics of the BQP model for the MSCP. 
Columns 1 to 3 give the instance names Instances along with the vertex num- 
ber V and edge number E of the graphs. Columns 4 and 5 show the number of 
colors K to be used and the number of variables N in the BQP formulation. 
Column 6 summarizes the best known results BKR from the previous liter- 
ature [61 IT1I1 111 15UJ51I2"5] . The columns under the heading of BQP-PR report 
our results of the BQP model solved by the PR algorithm: the best objective 
values Best, the average objective values Avr, the standard deviation a, the 
average time T\,_ avr (in seconds) to reach the best objective value Best over 
20 runs, and the average time (in seconds) Tavr consumed to reach the best 



http://mat.gsia.cmu.edu/COLOR/instances.html 
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Table 1 Computational statistics of the 


BQP-PR approach for the MSCP 






Instances 


V 


E 


K 


N 


BKR 




BQP-PR 








Best 


Avr 


a 


n_ avr 


Tav R 


mycicl3 


11 


20 


6 


66 


21* 


21 


21.0 


0.0 


< i 


< 1 


mycicl4 


23 


71 


7 


161 


45* 


45 


45.0 


0.0 


< i 


< 1 


mycicl5 


47 


23(5 


8 


376 


93* 


93 


93.0 


0.0 


2 


2 


mycicl6 


95 


755 


10 


950 


189* 


189 


189.0 


0.0 


497 


407 


mycicl7 


191 


2360 


10 


1910 


381 


381 


384.9 


3.7 


70 


384 


anna 


138 


493 


13 


1794 


276* 


276 


276.3 


0.4 


870 


721 


david 


87 


406 


13 


1131 


237* 


237 


237.0 


0.0 


524 


524 


huck 


74 


301 


13 


962 


243* 


243 


243.0 


0.0 


3 


3 


jean 


80 


254 


12 


960 


217* 


217 


217.0 


0.0 


79 


79 


quccn5.5 


25 


160 


7 


175 


75* 


75 


75.0 


0.0 


< 1 


< 1 


quccn6.6 


36 


290 


9 


324 


138* 


138 


138.0 


0.0 


6 





quccn7.7 


49 


476 


9 


441 


196* 


196 


196.0 


0.0 


6 


6 


quccn8.8 


64 


728 


10 


640 


291* 


291 


298.1 


5.1 


1064 


780 


gamcsl20 


120 


638 


10 


1200 


443* 


443 


446.5 


3.2 


755 


880 


milos250 


128 


387 


10 


1280 


325* 


325 


328.6 


2.3 


777 


1704 


milcs500 


128 


1170 


22 


2816 


705* 


713 


722.5 


6.4 


653 


1942 


DSJC125.1 


125 


736 


8 


1000 


326 


329 


338.5 


6.3 


1684 


7115 


DSJC125.5 


125 


3891 


22 


2750 


1015 


1050 


1082.6 


20.2 


32924 


15186 


DSJC125.9 


125 


6961 


50 


6250 


2511 


2529 


2573.8 


26.2 


34801 


24650 


DSJC250.1 


250 


3218 


12 


3000 


977 


1027 


1062.9 


16.8 


8893 


17206 


DSJC250.5 


250 


15668 


35 


8750 


3246 


3604 


3724.9 


59.1 


27009 


26065 


DSJC250.9 


250 


27897 


80 


20000 


8286 


8604 


8869.4 


122.2 


70737 


65673 


DSJC500.1 


500 


12458 


16 


8000 


2850 


3152 


3234.1 


41.7 


42447 


59241 


Average 












0.0160 


0.0283 


13.63 


9730.5 


9681.1 



164 objective value obtained in each run. The last row shows the average perfor- 
ms mance over a total of 23 tested instances in terms of the deviations of the best 
lee and average solution values from the BKR. Notice that the results marked in 

167 bold in the Best column indicate that BQP-PR reaches the BKR on these 

168 instances. 

mo We also applied CPLEX V12.2 to the linear model ((T|) fScctionl2~Tj) to solve 

™ these MSCP instances. CPLEX was successful in finding an optimal solution 

171 for 15 instances (marked with an asterisk), but terminated abnormally for the 

172 remaining 8 instances due to excess requirements of memory. 

173 From Table [TJ we observe that our BQP-PR approach is able to reach 

174 the best known results for 15 out of 23 instances, among which 14 results 

175 are known to be optimal values. The table discloses that only the few best 

176 algorithms that are specifically tailored to the MSCP can compete with this 

177 performance (see also Section |4~31 Table [2]). Moreover, the best and average 

178 solution values obtained by BQP-PR are very close to the best known results, 

179 with average deviations of 0.0160 and 0.0283, respectively over the set of the 
i8o benchmark instances. 

i8i A further observation is that our BQP-PR approach is quite robust to reach 

182 optimal or best known solution values within a short period of time for the 

183 instances with N < 2000 variables in comparison with a slow convergence for 

184 instances with many more BQP variables. We will further discuss this point 

185 in Section |4~H 
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Table 2 Comparison between 


the BQP-PR approach and other specific MSCP a] 


gorithms 


Instances 


BKR 


BQP-PR 


HLS 7 


MRLF 15 


PGAI18I 


TS6 


EXSCOL 25 


MDS(5) 11 


mycicl3 


21* 


21 


21 


21 


21 


- 


21 


21 


mycicl4 


45* 


45 


45 


45 


45 


- 


45 


45 


myciel5 


93* 


93 


93 


93 


93 


- 


93 


93 


mycicl6 


189* 


189 


189 


189 


189 


- 


189 


189 


myciel7 


381 


381 


381 


381 


382 


- 


381 


381 


anna 


276* 


276 


- 


277 


281 


- 


283 


276 


david 


237* 


237 


- 


241 


243 


- 


237 


237 


huck 


243* 


243 


243 


244 


243 


- 


243 


243 


jean 


217* 


217 


- 


217 


218 


- 


217 


217 


quccn5.5 


75* 


75 


- 


75 


75 


- 


75 


75 


quccn6.6 


138* 


138 


138 


138 


138 


- 


150 


138 


queen7.7 


196* 


196 


- 


196 


196 


- 


196 


196 


quccn8.8 


291* 


291 


- 


303 


302 


- 


291 


291 


gamcsl20 


443* 


443 


446 


446 


460 


- 


443 


443 


miles250 


325* 


325 


343 


334 


347 


- 


328 


325 


milcs500 


705* 


713 


755 


715 


762 


- 


709 


712 


DSJC125.1 


326 


329 


- 


352 


- 


344 


326 


326 


DSJC125.5 


1015 


1050 


- 


1141 


- 


1103 


1017 


1015 


DSJC125.9 


2511 


2529 


- 


2653 


- 


2631 


2512 


2511 


DSJC250.1 


977 


1027 


- 


1068 


- 


1046 


985 


977 


DS.IC250.5 


3246 


3604 


- 


3658 


- 


3779 


3246 


3281 


DSJC250.9 


8286 


8604 


- 


8942 


- 


9198 


8286 


8412 


DSJC500.1 


2850 


3152 


- 


3229 


- 


3205 


2850 


2951 



4.3 Comparison with other special purpose algorithms for the MSCP 
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In order to further evaluate our BQP-PR approach, we show a comparison of 
the proposed approach with several special purpose algorithms for the MSCP. 
These algorithms include a hybrid local search algorithm HLS [7] , an advanced 
recursive largest first algorithm MRLF [TS], a parallel genetic algorithm PGA 
[15] . a tabu search algorithm TS [§], a very recent independent set extrac- 
tion based heuristic EXSCOL [25 and a recent local search heuristic MSD(5) 
[11) . Given that the reference algorithms either utilize different termination 
conditions or do not report the computing time, we base the comparison on 
solution quality. Table [H presents the best objective values obtained by each 
algorithm (BQP-PR, HLS, MRLF, PGA, TS, EXSCOL and MSD(5), respec- 
tively) where the best solution values among them are marked in bold. The 
results for HLS, PGA and TS which are unavailable are marked with "-". 
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As we can observe in Tabled the proposed BQP-PR approach outperforms 
HLS, MRLF, PGA and TS in terms of the best solution values. Specifically, 
BQP-PR finds better solutions than HLS, MRLF, PGA and TS for 3, 14, 8 
and 7 instances, respectively. Finally, BQP-PR performs less well compared 
with the most effective MSCP heuristics EXSCOL and MDS(5). From this 
experiment, we conclude that our BQP approach combined with the PR algo- 
rithm for tacking the MSCP constitutes an interesting alternative to specific 
algorithms tailored to this problem. 
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4.4 Discussion 



In this section, we discuss some limitations of the proposed BQP-PR approach 
for solving the MSCP. First, the proposed method may require considerable 
computing time to reach its best solutions for large graph instances (see column 
Tavr hi Table [1]). This can be partially explained by the fact that the number 
of the BQP variables (equaling V ■ K where V is the number of vertices and 
K is the number of colors) sharply increases with the growth of V and K. 

214 Additionally, at present our approach is not able to solve graph instances with 

215 BQP variables surpassing the threshold value of 20,000 because of the memory 

216 limitation. These obstacles could be overcome by designing more effective data 
21? structures used by the BQP algorithms. 



5 Conclusion 

We have investigated the possibility of solving the NP-hard minimum sum col- 
oring problem (MSCP) via binary quadratic programming (BQP). We have 
shown how the MSCP can be recast into the BQP model and explained the 
key ideas of the Path-Relinking algorithm designed for the general BQP. Ex- 
periments on a set of benchmark instances demonstrate that this general BQP 
approach is able to reach competitive solutions compared with several special 
purpose MSCP algorithms even though considerable computing time may be 
required. 
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